package fem2.element.crack_sda;

import inf.jlinalg.Array1DMatrix;
import inf.jlinalg.ArrayVector;
import inf.jlinalg.IVector;
import inf.math.LinAlg;
import math2.MathUtilities;
import fem2.Debugger;
import fem2.MaterialModel;
import fem2.MeshPart;
import fem2.Node;
import fem2.enu.EchoLevelType;
import fem2.enu.StressType;
import fem2.material.cohesive_law.CohesiveLaw;

/**
 * CST element using Symmetric formulation, make the angle between sigma1 and
 * gradN as small as possible
 * 
 * @author hbui
 * 
 */
public class SOSConstantJumpCSTElement2 extends SOSConstantJumpCSTElement {

	public SOSConstantJumpCSTElement2(MeshPart mp, MaterialModel mm, CohesiveLaw cl) {
		super(mp, mm, cl);
	}

	public SOSConstantJumpCSTElement2(MeshPart mp, MaterialModel mm, CohesiveLaw cl,
			StressType stressType) {
		super(mp, mm, cl, stressType);
	}

	protected int posNode;

	@Override
	protected boolean checkOnsetCriteria() {
		IVector sigma = getStress(0, stressType);

		double oxx = sigma.get(0);
		double oyy = sigma.get(1);
		double oxy = sigma.get(2);
		double t1 = (oxx + oyy) / 2;
		double t2 = Math.sqrt(Math.pow((oxx - oyy) / 2, 2) + Math.pow(oxy, 2));
		double o1 = t1 + t2;

		if (o1 < getCohesiveLaw().getTensileStrength()) {
			return false;
		} else {
			MeshPart mp = getMeshPart();
			int n = mp.getBasis().getDimension();
			int dim = mp.getDimension();
			int ntrs = dim * (dim + 1) / 2;
			alpha = new ArrayVector(dim);

			/*
			 * realize crack normal
			 */
			double n1, n2;
			if (Math.abs(oxy) < MathUtilities.tol14) {
				if (oxx > oyy) {
					n1 = 1;
					n2 = 0;
				} else {
					n1 = 0;
					n2 = 1;
				}
			} else {
				double g = Math.sqrt(Math.pow(oxx - o1, 2) + Math.pow(oxy, 2));
				n1 = -oxy / g;
				n2 = (oxx - o1) / g;
			}

			/*
			 * compute B operator
			 */
			Array1DMatrix B = Array1DMatrix.createColumnMajor(ntrs, dim * n);
			double[] xi = getIntegrationPoint(0);
			computeB(mp, B, xi);

			/*
			 * search the gradient which make smallest angle with crack normal
			 */
			double cos = 0.0;
			posNode = 0;
			for (int i = 0; i < 3; i++) {
				double gradx = B.get(0, i * dim);
				double grady = B.get(1, 1 + i * dim);
				double normgrad = Math.sqrt(gradx * gradx + grady * grady);
				gradx /= normgrad;
				grady /= normgrad;
				double tmp = Math.abs(n1 * gradx + n2 * grady);
				if (tmp > cos) {
					cos = tmp;
					posNode = i;
					/*
					 * choose normal to make smallest angle with gradient vector
					 */
					double tmp1 = n1 * gradx + n2 * grady;
					normal[0] = n1 * Math.signum(tmp1);
					normal[1] = n2 * Math.signum(tmp1);
					/*
					 * compute crack length
					 */
					crackMeasure = elementMeasure * normgrad;
				}
			}
			/*
			 * compute tangent
			 */
			tangent1[0] = -normal[1];
			tangent1[1] = normal[0];

			/*
			 * compute crack length and two intersections
			 */
			int[] negNodes;
			if (posNode == 0) {
				negNodes = new int[] { 1, 2 };
			} else if (posNode == 1) {
				negNodes = new int[] { 0, 2 };
			} else if (posNode == 2) {
				negNodes = new int[] { 0, 1 };
			} else {
				throw new Error("positive node is wrong");
			}
			/*
			 * realize two edges intersect with crack
			 */
			int[] edge1 = new int[] { posNode, negNodes[0] };
			int[] edge2 = new int[] { posNode, negNodes[1] };

			/*
			 * set the crack center as element center and compute intersection
			 * with edge 1
			 */
			Node[] nodes = mp.getNodes();
			double[] N1 = nodes[edge1[0]].getPosition();
			double[] N2 = nodes[edge1[1]].getPosition();
			double xi1 = N1[0], yi1 = N1[1];
			double xi2 = N2[0], yi2 = N2[1];
			double a1 = yi2 - yi1;
			double b1 = xi1 - xi2;
			double c1 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
			double[] cen = mp.computeCenter();
			double a = normal[0];
			double b = normal[1];
			double c = -(a * cen[0] + b * cen[1]);
			double[] X1 = MathUtilities.computeIntersection(a, b, c, a1, b1, c1);
			/*
			 * compute intersection with edge 2
			 */
			N1 = nodes[edge2[0]].getPosition();
			N2 = nodes[edge2[1]].getPosition();
			xi1 = N1[0];
			yi1 = N1[1];
			xi2 = N2[0];
			yi2 = N2[1];
			double a2 = yi2 - yi1;
			double b2 = xi1 - xi2;
			double c2 = yi1 * (xi2 - xi1) - xi1 * (yi2 - yi1);
			double[] X2 = MathUtilities.computeIntersection(a, b, c, a2, b2, c2);

			/*
			 * compute temporary length
			 */
			double tmpLength = LinAlg.diffNorm2(2, X1, X2);
			/*
			 * realize actual crack center by scaling
			 */
			double scale = crackMeasure / tmpLength;
			double[] crackCenter = new double[2];
			double[] posPos = nodes[posNode].getPosition();
			crackCenter[0] = (cen[0] - posPos[0]) * scale + posPos[0];
			crackCenter[1] = (cen[1] - posPos[1]) * scale + posPos[1];

			/*
			 * compute crack boundary
			 */
			c = -(a * crackCenter[0] + b * crackCenter[1]);
			X1 = MathUtilities.computeIntersection(a, b, c, a1, b1, c1);
			X2 = MathUtilities.computeIntersection(a, b, c, a2, b2, c2);
			crackBoundary.add(X1);
			crackBoundary.add(X2);

			/*
			 * debugging part
			 */
			if (getEchoLevel().contain(EchoLevelType.AUX4)) {
				Debugger.watch("crack length = ", crackMeasure);
				Debugger.watch("crack normal = ", normal);
				Debugger.watch("element area = ", elementMeasure);
			}

			return true;
		}
	}
}
